#include "spline_lib.h"

double bvalue ( double *t, double *bcoef, int n, int k, double x, int jderiv )
{
/*  from  * a practical guide to splines *  by c. de boor    
calls  interv
c
calculates value at  x  of  jderiv-th derivative of spline from b-repr.
c  the spline is taken to be continuous from the right, EXCEPT at the
c  rightmost knot, where it is taken to be continuous from the left.
c
c******  i n p u t ******
c  t, bcoef, n, k......forms the b-representation of the spline  f  to
c        be evaluated. specifically,
c  t.....knot sequence, of length  n+k, assumed nondecreasing.
c  bcoef.....b-coefficient sequence, of length  n .
c  n.....length of  bcoef  and dimension of spline(k,t),
c        a s s u m e d  positive .
c  k.....order of the spline .
c
c  w a r n i n g . . .   the restriction  k .le. kmax (=20)  is imposed
c        arbitrarily by the dimension statement for  aj, dl, dr  below,
c        but is  n o w h e r e  c h e c k e d  for.
c
c  x.....the point at which to evaluate .
c  jderiv.....integer giving the order of the derivative to be evaluated
c        a s s u m e d  to be zero or positive.
c
c******  o u t p u t  ******
c  bvalue.....the value of the (jderiv)-th derivative of  f  at  x .
c
c******  m e t h o d  ******
c     The nontrivial knot interval  (t[i),t[i+1])  containing  x  is lo-
c  cated with the aid of  interv . The  k  b-coeffs of  f  relevant for
c  this interval are then obtained from  bcoef (or taken to be zero if
c  not explicitly available) and are then differenced  jderiv  times to
c  obtain the b-coeffs of  (d**jderiv)f  relevant for that interval.
c  Precisely, with  j = jderiv, we have from x.(12) of the text that
c
c     (d**j)f  =  sum ( bcoef[.,j)*b(.,k-j,t) )
c
c  where
c                   / bcoef[.),                     ,  j .eq. 0
c                   /
c    bcoef[.,j)  =  / bcoef[.,j-1] - bcoef[.-1,j-1]
c                   / ----------------------------- ,  j .gt. 0
c                   /    (t[.+k-j] - t[.))/(k-j)
c
c     Then, we use repeatedly the fact that
c
c    sum ( a(.)*b(.,m,t)(x) )  =  sum ( a(.,x)*b(.,m-1,t)(x) )
c  with
c                 (x - t[.))*a(.) + (t[.+m-1] - x)*a(.-1]
c    a(.,x)  =    ---------------------------------------
c                 (x - t[.))      + (t[.+m-1] - x)
c
c  to write  (d**j)f(x)  eventually as a linear combination of b-splines
c  of order  1 , and the coefficient for  b(i,1,t)(x)  must then be the
c  desired number  (d**j)f(x). (see x.(17)-(19) of text).
c 
*/
      const int kmax = 20;
      int   i,ilo,imk,j,jc,jcmin,jcmax,jj,kmj,km1,mflag
                    ,nmi,jdrvp1;

      double   aj[21],dl[21],dr[21], fkmj, bvalue_ret;

      bvalue_ret = 0.;
      if (jderiv >= k)                goto label_99;
/*
c
c  *** Find  i   s.t.   1 .le. i .lt. n+k   and   t[i) .lt. t[i+1]   and
c      t[i) .le. x .lt. t[i+1] . If no such i can be found,  x  lies
c      outside the support of  the spline  f , hence  bvalue = 0.
c      (The asymmetry in this choice of  i  makes  f  rightcontinuous, except
c      at  t[n+k] where it is leftcontinuous.) */
      interv ( t, n+k, x, &i, &mflag );
      if (mflag != 0)                 goto label_99;
/*
c  *** if k = 1 (and jderiv = 0), bvalue = bcoef[i].
*/
      km1 = k - 1;
      if (km1 > 0)                   goto label_1;
      bvalue_ret = bcoef[i];
                                        goto label_99;
/*
c
c  *** store the k b-spline coefficients relevant for the knot interval
c     (t[i),t[i+1]) in aj[1],...,aj[k) and compute dl[j] = x - t[i+1-j],
c     dr[j] = t[i+j] - x, j=1,...,k-1 . set any of the aj not obtainable
c     from input to zero. set any t.s not obtainable equal to t[1] or
c     to t[n+k) appropriately.
*/
label_1:
	  jcmin = 1;
      imk = i - k;
      if (imk >= 0)                   goto label_8;
      jcmin = 1 - imk;
	  for(j=1; j<=i; j++)
		  dl[j] = x - t[i+1-j];
      for(j=i; j<= km1; j++)
          aj[k-j] = 0., dl[j] = dl[i];
    
      
                                        goto label_10;
label_8:
	  for(j=1; j<=km1; j++)
		dl[j] = x - t[i+1-j];
	
label_10: 
	  jcmax = k;
      nmi = n - i;
      if (nmi >= 0)                   goto label_18;
      jcmax = k + nmi;
	  for(j=1; j<= jcmax; j++)
		  dr[j] = t[i+j] - x;

      for(j=jcmax; j<=km1; j++)
		  aj[j+1] = 0.,
	      dr[j] = dr[jcmax];

      
                                        goto label_20;
   
label_18:	  
		for(j=1; j<= km1; j++)
		  dr[j] = t[i+j] - x;								
										
   
label_20:    
		for(jc=jcmin; jc<=jcmax; jc++)
           aj[jc] = bcoef[imk + jc];
		  
		   
/*		   
c
c               *** difference the coefficients  jderiv  times.
*/
      if (jderiv == 0)                goto label_30;
	  for(j=1; j<=jderiv; j++){
		  kmj = k-j;
		  fkmj = (double) kmj;
		  ilo = kmj;
		  for(jj=1; jj<=kmj; jj++){
			  aj[jj] = ((aj[jj+1] - aj[jj])/(dl[ilo] + dr[jj]))*fkmj;
              ilo = ilo - 1;
		  }
	  }
/*
c
c  *** compute value at  x  in (t[i),t[i+1]) of jderiv-th derivative,
c     given its relevant b-spline coeffs in aj[1],...,aj[k-jderiv).
*/
label_30:  
    if (jderiv == km1)              goto label_39;
      jdrvp1 = jderiv + 1;
    for(j=jdrvp1; j<=km1; j++){
		kmj = k-j;
		ilo = kmj;
		for(jj=1; jj<=kmj; jj++){
			aj[jj] = (aj[jj+1]*dl[ilo] + aj[jj]*dr[jj])/(dl[ilo]+dr[jj]);
			ilo = ilo - 1;
		}
	}
label_39:
	bvalue_ret = aj[1];
label_99:
	return bvalue_ret;
}



void interv ( double * xt, int lxt, double x, int * left, int * mflag )
{
/*
c  from  * a practical guide to splines *  by C. de Boor    
computes  left = max( i :  xt(i) .lt. xt(lxt) .and.  xt(i) .le. x )  .
c
c******  i n p u t  ******
c  xt.....a real sequence, of length  lxt , assumed to be nondecreasing
c  lxt.....number of terms in the sequence  xt .
c  x.....the point whose location with respect to the sequence  xt  is
c        to be determined.
c
c******  o u t p u t  ******
c  left, mflag.....both integers, whose value is
c
c   1     -1      if               x .lt.  xt(1)
c   i      0      if   xt(i)  .le. x .lt. xt(i+1)
c   i      0      if   xt(i)  .lt. x .eq. xt(i+1) .eq. xt(lxt)
c   i      1      if   xt(i)  .lt.        xt(i+1) .eq. xt(lxt) .lt. x
c
c        In particular,  mflag = 0  is the 'usual' case.  mflag .ne. 0
c        indicates that  x  lies outside the CLOSED interval
c        xt(1) .le. y .le. xt(lxt) . The asymmetric treatment of the
c        intervals is due to the decision to make all pp functions cont-
c        inuous from the right, but, by returning  mflag = 0  even if
C        x = xt(lxt), there is the option of having the computed pp function
c        continuous from the left at  xt(lxt) .
c
c******  m e t h o d  ******
c  The program is designed to be efficient in the common situation that
c  it is called repeatedly, with  x  taken from an increasing or decrea-
c  sing sequence. This will happen, e.g., when a pp function is to be
c  graphed. The first guess for  left  is therefore taken to be the val-
c  ue returned at the previous call and stored in the  l o c a l  varia-
c  ble  ilo . A first check ascertains that  ilo .lt. lxt (this is nec-
c  essary since the present call may have nothing to do with the previ-
c  ous call). Then, if  xt(ilo) .le. x .lt. xt(ilo+1), we set  left =
c  ilo  and are done after just three comparisons.
c     Otherwise, we repeatedly double the difference  istep = ihi - ilo
c  while also moving  ilo  and  ihi  in the direction of  x , until
c                      xt(ilo) .le. x .lt. xt(ihi) ,
c  after which we use bisection to get, in addition, ilo+1 = ihi .
c  left = ilo  is then returned.
c
*/
      int  ihi, istep, middle;

      static int  ilo; 
      
      if(ilo <1 || ilo>lxt) ilo = 1;
      ihi = ilo + 1;
      if (ihi < lxt)                  goto label_20;
         if (x >= xt[lxt])            goto label_110;
         if (lxt <= 1)                goto label_90;
         ilo = lxt - 1;
         ihi = lxt;

label_20:
      if (x >= xt[ihi])               goto label_40;
      if (x >= xt[ilo])               goto label_100;
/*
              **** now x .lt. xt(ilo) . decrease  ilo  to capture  x .
*/
      istep = 1;
label_31:
         ihi = ilo;
         ilo = ihi - istep;
         if (ilo <= 1)                  goto label_35;
         if (x >= xt[ilo])              goto label_50;
         istep = istep*2;
										goto label_31;
label_35:
	  ilo = 1;
      if (x < xt[1])                    goto label_90;
                                        goto label_50;
/*
              **** now x .ge. xt(ihi) . increase  ihi  to capture  x .
*/
label_40:   
   istep = 1;
label_41:
         ilo = ihi;
         ihi = ilo + istep;
         if (ihi >= lxt)              goto label_45;
         if (x < xt[ihi])             goto label_50;
         istep = istep*2;
                                      goto label_41;
label_45:
      if (x >= xt[lxt])               goto label_110;
      ihi = lxt;
/*
c           **** now xt(ilo) .le. x .lt. xt(ihi) . narrow the interval.
*/
label_50:
      middle = (ilo + ihi)/2;
      if (middle == ilo)              goto label_100;
/*
c     note. it is assumed that middle = ilo in case ihi = ilo+1 .
*/
      if (x < xt[middle])              goto label_53;
         ilo = middle;
                                       goto label_50;
label_53:
      ihi = middle;
                                       goto label_50;
/* c**** set output and return.  */
label_90:
      *mflag = -1;
      *left = 1;
                                        return;
label_100:  
	  *mflag = 0;
      *left = ilo;
                                        return;
label_110:  
	  *mflag = 1;
      if (x == xt[lxt]) *mflag = 0;
      *left = lxt;
label_111:
      if (*left == 1)                  return;
      *left = *left - 1;
      if (xt[*left] < xt[lxt])        return;
          goto label_111;
}



void bsplvb (double * t, int jhigh, int index, double x, int left, double * biatx )
{
/*
c  from  * a practical guide to splines *  by c. de boor    
calculates the value of all possibly nonzero b-splines at  x  of order
c
c               jout  =  max( jhigh , (j+1)*(index-1) )
c
c  with knot sequence  t .
c
c******  i n p u t  ******
c  t.....knot sequence, of length  left + jout  , assumed to be nonde-
c        creasing.  a s s u m p t i o n . . . .
c                       t(left)  .lt.  t(left + 1)   .
c   d i v i s i o n  b y  z e r o  will result if  t(left) = t(left+1)
c  jhigh,
c  index.....integers which determine the order  jout = max(jhigh,
c        (j+1)*(index-1))  of the b-splines whose values at  x  are to
c        be returned.  index  is used to avoid recalculations when seve-
c        ral columns of the triangular array of b-spline values are nee-
c        ded (e.g., in  bsplpp  or in  bsplvd ). precisely,
c                     if  index = 1 ,
c        the calculation starts from scratch and the entire triangular
c        array of b-spline values of orders 1,2,...,jhigh  is generated
c        order by order , i.e., column by column .
c                     if  index = 2 ,
c        only the b-spline values of order  j+1, j+2, ..., jout  are ge-
c        nerated, the assumption being that  biatx , j , deltal , deltar
c        are, on entry, as they were on exit at the previous call.
c           in particular, if  jhigh = 0, then  jout = j+1, i.e., just
c        the next column of b-spline values is generated.
c
c  w a r n i n g . . .  the restriction   jout .le. jmax (= 20)  is im-
c        posed arbitrarily by the dimension statement for  deltal  and
c        deltar  below, but is  n o w h e r e  c h e c k e d  for .
c
c  x.....the point at which the b-splines are to be evaluated.
c  left.....an integer chosen (usually) so that
c                  t(left) .le. x .le. t(left+1)  .
c
c******  o u t p u t  ******
c  biatx.....array of length  jout , with  biatx(i)  containing the val-
c        ue at  x  of the polynomial of order  jout  which agrees with
c        the b-spline  b(left-jout+i,jout,t)  on the interval (t(left),
c        t(left+1)) .
c
c******  m e t h o d  ******
c  the recurrence relation
c
c                       x - t(i)              t(i+j+1) - x
c     b(i,j+1)(x)  =  -----------b(i,j)(x) + ---------------b(i+1,j)(x)
c                     t(i+j)-t(i)            t(i+j+1)-t(i+1)
c
c  is used (repeatedly) to generate the (j+1)-vector  b(left-j,j+1)(x),
c  ...,b(left,j+1)(x)  from the j-vector  b(left-j+1,j)(x),...,
c  b(left,j)(x), storing the new values in  biatx  over the old. the
c  facts that
c            b(i,1) = 1  if  t(i) .le. x .lt. t(i+1)
c  and that
c            b(i,j)(x) = 0  unless  t(i) .le. x .lt. t(i+j)
c  are used. the particular organization of the calculations follows al-
c  gorithm  (8)  in chapter x of the text.
c
      
current fortran standard makes it impossible to specify the length of
c  t  and of  biatx  precisely without the introduction of otherwise
c  superfluous additional arguments.
*/

      static int  j;
	  static double deltal[21],deltar[21]; 
	  const int jmax = 20;
      int   i, jp1;
      double  saved, term;          
      
	  if(j <= 0) j=1;	
      if(index == 1) 
		  goto label_10;
	  else
		  goto label_20;

label_10:                          
      j = 1;
      biatx[1] = 1.;
      if (j >= jhigh)                 goto label_99;

label_20:
         jp1 = j + 1;
         deltar[j] = t[left+j] - x;
         deltal[j] = x - t[left+1-j];
         saved = 0.;
		 for(i=1; i<=j; i++){
            term = biatx[i]/(deltar[i] + deltal[jp1-i]);
            biatx[i] = saved + deltar[i]*term;
            saved = deltal[jp1-i]*term;
		 }
         biatx[jp1] = saved;
         j = jp1;
         if (j < jhigh)              goto label_20;
label_99:
      return;
}
    

void bsplvd (double * t, int k, double x, int left, 
			 double * dbiatx, int nderiv )
{
/*
c  from  * a practical guide to splines *  by c. de Boor (7 may 92)    
calls bsplvb
calculates value and deriv.s of all b-splines which do not vanish at x
c
c******  i n p u t  ******
c  t     the knot array, of length left+k (at least)
c  k     the order of the b-splines to be evaluated
c  x     the point at which these values are sought
c  left  an integer indicating the left endpoint of the interval of
c        interest. the  k  b-splines whose support contains the interval
c               (t(left), t(left+1))
c        are to be considered.
c  a s s u m p t i o n  - - -  it is assumed that
c               t(left) .lt. t(left+1)
c        division by zero will result otherwise (in  b s p l v b ).
c        also, the output is as advertised only if
c               t(left) .le. x .le. t(left+1) .
c  nderiv   an integer indicating that values of b-splines and their
c        derivatives up to but not including the  nderiv-th  are asked
c        for. ( nderiv  is replaced internally by the integer  m h i g h
c        in  (1,k)  closest to it.)
c
c******  w o r k   a r e a  ******
c  a     an array of order (k,k), to contain b-coeff.s of the derivat-
c        ives of a certain order of the  k  b-splines of interest.
c
c******  o u t p u t  ******
c  dbiatx   an array of order (k,nderiv). its entry  (i,m)  contains
c        value of  (m-1)st  derivative of  (left-k+i)-th  b-spline of
c        order  k  for knot sequence  t , i=1,...,k, m=1,...,nderiv.
c
c******  m e t h o d  ******
c  values at  x  of all the relevant b-splines of order k,k-1,...,
c  k+1-nderiv  are generated via  bsplvb  and stored temporarily in
c  dbiatx .  then, the b-coeffs of the required derivatives of the b-
c  splines of interest are generated by differencing, each from the pre-
c  ceding one of lower order, and combined with the values of b-splines
c  of corresponding order in  dbiatx  to produce the desired values .
c
*/
      int i, ideriv, il, j, jlow, jp1mid, kp1, kp1mm
                            ,ldummy, m, mhigh;
      double  *a, factor, fkp1mm, sum, *temp ;
                          
      a = (double *) calloc( (k+1)*(k+1), sizeof(double));
	  temp = (double *) calloc( (k+1), sizeof(double));
	  if(nderiv > k) mhigh =k;
	  else mhigh = nderiv;
/*    mhigh is usually equal to nderiv. */
      kp1 = k+1;
      bsplvb(t,kp1-mhigh,1,x,left, temp);
	  for(i=1; i<=k; i++)
		  dbiatx[i+k*1] = temp[i];
      if (mhigh == 1)                 goto label_99;
/*
c     the first column of  dbiatx  always contains the b-spline values
c     for the current order. these are stored in column k+1-current
c     order  before  bsplvb  is called to put values for the next
c     higher order on top of it.
*/
      ideriv = mhigh;
	  for(m=2; m<=mhigh; m++){
      
         jp1mid = 1;
		 for(j=ideriv; j<=k; j++){
            dbiatx[j+k*ideriv] = dbiatx[jp1mid+k*1];
            jp1mid = jp1mid + 1;
		 }
         ideriv = ideriv - 1;
         bsplvb(t,kp1-ideriv,2,x,left, temp);
         for(i=1; i<=k; i++)
		  dbiatx[i+k*1] = temp[i];
	  }
/*
c     at this point,  b(left-k+i, k+1-j)(x) is in  dbiatx(i,j) for
c     i=j,...,k and j=1,...,mhigh ('=' nderiv). in particular, the
c     first column of  dbiatx  is already in final form. to obtain cor-
c     responding derivatives of b-splines in subsequent columns, gene-
c     rate their b-repr. by differencing, then evaluate at  x.
*/
      jlow = 1;
	  for(i=1; i<=k; i++){
		for(j=jlow; j<=k; j++)
          a[j+i*k] = 0.;
      
        jlow = i;
        a[i+k*i] = 1.;
	  }
/*
c     at this point, a(.,j) contains the b-coeffs for the j-th of the
c     k  b-splines of interest here.
*/
     
	for(m=2; m<=mhigh; m++){
         kp1mm = kp1 - m;
         fkp1mm = (double)kp1mm;
         il = left;
         i = k;
/*
c        for j=1,...,k, construct b-coeffs of  (m-1)st  derivative of
c        b-splines from those for preceding derivative by differencing
c        and store again in  a(.,j) . the fact that  a(i,j) = 0  for
c        i .lt. j  is used.
*/
        
		for(ldummy=1; ldummy<= kp1mm; ldummy++){
            factor = fkp1mm/(t[il+kp1mm] - t[il]);
/*
c           the assumption that t(left).lt.t(left+1) makes denominator
c           in  factor  nonzero.
*/
            
            for(j=1; j<=i; j++)
               a[i+k*j] = (a[i+k*j] - a[i-1+k*j])*factor;
            il = il - 1;
            i = i - 1;
		}
/*
c        for i=1,...,k, combine b-coeffs a(.,i) with b-spline values
c        stored in dbiatx(.,m) to get value of  (m-1)st  derivative of
c        i-th b-spline (of interest here) at  x , and store in
c        dbiatx(i,m). storage of this value over the value of a b-spline
c        of order m there is safe since the remaining b-spline derivat-
c        ives of the same order do not use this value due to the fact
c        that  a(j,i) = 0  for j .lt. i .
*/
         
		 for(i=1; i<=k; i++){
            sum = 0.;
			jlow = m;
			if(i>m) jlow =i;
          
			for(j = jlow; j<=k; j++)
				sum = a[j+k*i]*dbiatx[j+k*m] + sum;
			dbiatx[i+k*m] = sum ;
		 }
	}
label_99:  
	free(a), free(temp);
   return;
}



struct SP_2D_TP * sp_2d_initialize(int order, int n1, int n2,
                                 int  * inform)
{
/*  Purpose: This function allocates and initializes all the pointers
    in the SP_2D_TP data structure.  

    Return Value:  A pointer to the newly initialized SP_2D_TP structure


    Arguments:

        order           Order of spline (input)
        n1              Domain dimension (input)
        n2              Range dimension (input)
        
        
        inform          Error code (output)


    Fatal Errors (which return a NULL pointer):

        inform          Error Condition

          1             order     <= 0
          2             n1        <= 0
          3             n2        <= 0
          4             allocation failure for some member of SP_2D_TP

    */

    struct SP_2D_TP * sp;
    int i;
    double * temp_ptr1;
    

    /* Initial error testing */
    *inform = 0;
    if(order <= 0)
        *inform = 1;
    else if(n1 <= 0)
        *inform = 2;
    else if(n2 <= 0)
        *inform = 3;
    if(*inform > 0) return NULL;



            /*  Allocate space for the SP_2D_TP structure  */
    sp = (struct SP_2D_TP *) CALLOC(1, sizeof(struct SP_2D_TP));
    if((void *)sp == NULL){
        *inform = 4;
        return NULL;
    }


            /*  Place numerical values in SP_2D_TP structure  */
    sp->order       = order;
    sp->n1		    = n1;
    sp->n2		    = n2;
    

	sp->coef = get_mat(n1, n2);
	if((void *)(sp->coef) == NULL){
        *inform = 4;
        return NULL;
    }

	sp->knot = (double **) CALLOC(2, sizeof(double *));
    if((void *)(sp->knot) == NULL){
        *inform = 4;
        return NULL;
    }
	
	 temp_ptr1   = (double *) CALLOC(n1+n2+2*order , sizeof(double));
    if((void *)temp_ptr1 == NULL){
        *inform = 4;
        return NULL;}

    for(i=0; i< 2; i++)
        (sp->knot)[i] = &(temp_ptr1[i*(n1+order)]);
	
	return sp;

}

struct SP_3D_TP * sp_3d_initialize(int order, int n1, int n2, int n3,
                                 int  * inform)
{
/*  Purpose: This function allocates and initializes all the pointers
    in the SP_2D_TP data structure.  

    Return Value:  A pointer to the newly initialized SP_2D_TP structure


    Arguments:

        order           Order of spline (input)
        n1              1st dimension (input)
        n2              2nd dimension (input)
		n3				3rd dimension (input)
        
        
        inform          Error code (output)


    Fatal Errors (which return a NULL pointer):

        inform          Error Condition

          1             order     <= 0
          2             n1        <= 0
          3             n2        <= 0
		  4				n3		  <= 0
          5             allocation failure for some member of SP_2D_TP

    */

    struct SP_3D_TP * sp;
    double * temp_ptr1;
    

    /* Initial error testing */
    *inform = 0;
    if(order <= 0)
        *inform = 1;
    else if(n1 <= 0)
        *inform = 2;
    else if(n2 <= 0)
        *inform = 3;
	else if(n3 <= 0)
		*inform = 4;
    if(*inform > 0) return NULL;



            /*  Allocate space for the SP_2D_TP structure  */
    sp = (struct SP_3D_TP *) CALLOC(1, sizeof(struct SP_3D_TP));
    if((void *)sp == NULL){
        *inform = 5;
        return NULL;
    }


            /*  Place numerical values in SP_2D_TP structure  */
    sp->order       = order;
    sp->n1		    = n1;
    sp->n2		    = n2;
	sp->n3			= n3;
    

	sp->coef = get_mat_3d(n1, n2, n3);
	if((void *)(sp->coef) == NULL){
        *inform = 5;
        return NULL;
    }

	sp->knot = (double **) CALLOC(3, sizeof(double *));
    if((void *)(sp->knot) == NULL){
        *inform = 5;
        return NULL;
    }
	
	 temp_ptr1   = (double *) CALLOC(n1+n2+n3+3*order , sizeof(double));
    if((void *)temp_ptr1 == NULL){
        *inform = 5;
        return NULL;}

    (sp->knot)[0] = &(temp_ptr1[0]);
	(sp->knot)[1] = &(temp_ptr1[n1+order]);
	(sp->knot)[2] = &(temp_ptr1[n1+n2+2*order]);
	
	return sp;

}

void sp_3d_initialize_without_calloc(struct SP_3D_TP *sp, int order, int n1, int n2, int n3,
                                 int  * inform)
{
/*  Purpose: This function allocates and initializes all the pointers
    in the SP_2D_TP data structure.  

    Return Value:  A pointer to the newly initialized SP_2D_TP structure


    Arguments:

        order           Order of spline (input)
        n1              1st dimension (input)
        n2              2nd dimension (input)
		n3				3rd dimension (input)
        
        
        inform          Error code (output)


    Fatal Errors (which return a NULL pointer):

        inform          Error Condition

          1             order     <= 0
          2             n1        <= 0
          3             n2        <= 0
		  4				n3		  <= 0
          5             allocation failure for some member of SP_2D_TP

    */

    double * temp_ptr1;
    

    /* Initial error testing */
    *inform = 0;
    if(order <= 0)
        *inform = 1;
    else if(n1 <= 0)
        *inform = 2;
    else if(n2 <= 0)
        *inform = 3;
	else if(n3 <= 0)
		*inform = 4;
    if(*inform > 0) return;



            /*  Place numerical values in SP_2D_TP structure  */
    sp->order       = order;
    sp->n1		    = n1;
    sp->n2		    = n2;
	sp->n3			= n3;
    

	sp->coef = get_mat_3d(n1, n2, n3);
	if((void *)(sp->coef) == NULL){
        *inform = 5;
        return;
    }

	sp->knot = (double **) CALLOC(3, sizeof(double *));
    if((void *)(sp->knot) == NULL){
        *inform = 5;
        return;
    }
	
	 temp_ptr1   = (double *) CALLOC(n1+n2+n3+3*order , sizeof(double));
    if((void *)temp_ptr1 == NULL){
        *inform = 5;
        return;}

    (sp->knot)[0] = &(temp_ptr1[0]);
	(sp->knot)[1] = &(temp_ptr1[n1+order]);
	(sp->knot)[2] = &(temp_ptr1[n1+n2+2*order]);
	
	return;

}


struct SP_1D * sp_1d_initialize(int order, int n,
                                 int  * inform)
{
/*  Purpose: This function allocates and initializes all the pointers
    in the SP_1D data structure.  

    Return Value:  A pointer to the newly initialized SP_1D structure.


    Arguments:

        order           Order of spline (input)
        n               Spline space dimension (input)
        
        
        inform          Error code (output)


    Fatal Errors (which return a NULL pointer):

        inform          Error Condition

          1             order  <= 0
          2             n      <= 0
          4             allocation failure for some member of SP_2D_TP

    */

    struct SP_1D * sp;
    double * temp_ptr1;
    

    /* Initial error testing */
    *inform = 0;
    if(order <= 0)
        *inform = 1;
    else if(n <= 0)
        *inform = 2;
    if(*inform > 0) return NULL;



            /*  Allocate space for the SP_2D_TP structure  */
    sp = (struct SP_1D *) CALLOC(1, sizeof(struct SP_1D));
    if((void *)sp == NULL){
        *inform = 4;
        return NULL;
    }


            /*  Place numerical values in SP_2D_TP structure  */
    sp->order       = order;
    sp->n		    = n;  

            /*  Allocate space for the coefficients
                and set the SP_1D pointers  */
    temp_ptr1   = (double *) CALLOC(n, sizeof(double));
    if((void *)temp_ptr1 == NULL){
        *inform = 4;
        return NULL;
    }

    
    
    (sp->coef)= temp_ptr1;

            /*  Allocate space for the knots
                and set the SP_2D_TP pointers  */
   
    temp_ptr1   = (double *) CALLOC(n+order , sizeof(double));
    if((void *)temp_ptr1 == NULL){
        *inform = 4;
        return NULL;}

	sp->knot = temp_ptr1;

    
	
	return sp;

}

void sp_1d_set_knots(struct SP_1D * sp, double a, double b)
{
/*  Purpose: This function returns equally spaced knots for the 1d 
             spline.  The x-knots partition the interval [a, b]
			 [a,b].

    Return Value:  void


    Arguments:

        sp              A 1d spline (input)
        a               Left end point for the x-knots (input)
        b               Right end point for the x-knots (input)
		c               Left end point for the y-knots (input)
        d               Right end point for the y-knots (input)               
        
        


    There are no Fatal Errors 

    */
int i;
double temp;

for(i=0; i<sp->order; i++){
	sp->knot[i] = a;
    sp->knot[sp->n+i] = b;
}

temp = b-a;
for(i=sp->order; i<sp->n; i++)
	sp->knot[i] = a + ((i-sp->order+1)*temp)/(sp->n - sp->order+1);
}


void sp_2d_set_knots(struct SP_2D_TP * sp, double a, double b, 
				  double c, double d)
{
/*  Purpose: This function returns equally spaced knots for the 2d 
             tensor product spline.  The x-knots partition the interval
			 [a,b] while the y-knots partition the interval [c,d].
			 It is assumed that a<b and c<d.

    Return Value:  void


    Arguments:

        sp              A 2d tensor product spline (input)
        a               Left end point for the x-knots (input)
        b               Right end point for the x-knots (input)
		c               Left end point for the y-knots (input)
        d               Right end point for the y-knots (input)               
        
        


    There are no Fatal Errors 

    */
int i;
double temp;

for(i=0; i<sp->order; i++){
	sp->knot[0][i] = a;
    sp->knot[0][sp->n1+i] = b;
	sp->knot[1][i] = c;
	sp->knot[1][sp->n2+i] = d;
}

temp = b-a;
for(i=sp->order; i<sp->n1; i++)
	sp->knot[0][i] = a + ((i-sp->order+1)*temp)/(sp->n1 - sp->order+1);

temp = d-c;
for(i=sp->order; i<sp->n2; i++)
	sp->knot[1][i] = c + ((i-sp->order+1)*temp)/(sp->n2 - sp->order+1);
}

void sp_3d_set_knots(struct SP_3D_TP * sp, double a, double b, 
				  double c, double d, double s, double t)
{
/*  Purpose: This function returns equally spaced knots for the 2d 
             tensor product spline.  The x-knots partition the interval
			 [a,b] while the y-knots partition the interval [c,d].
			 It is assumed that a<b and c<d.

    Return Value:  void


    Arguments:

        sp              A 2d tensor product spline (input)
        a               Left end point for the x-knots (input)
        b               Right end point for the x-knots (input)
		c               Left end point for the y-knots (input)
        d               Right end point for the y-knots (input)               
        
        


    There are no Fatal Errors 

    */
int i;
double temp;

for(i=0; i<sp->order; i++){
	sp->knot[0][i] = a;
    sp->knot[0][sp->n1+i] = b;
	sp->knot[1][i] = c;
	sp->knot[1][sp->n2+i] = d;
	sp->knot[2][i] = s;
	sp->knot[2][sp->n3+i] = t;
}

temp = b-a;
for(i=sp->order; i<sp->n1; i++)
	sp->knot[0][i] = a + ((i-sp->order+1)*temp)/(sp->n1 - sp->order+1);

temp = d-c;
for(i=sp->order; i<sp->n2; i++)
	sp->knot[1][i] = c + ((i-sp->order+1)*temp)/(sp->n2 - sp->order+1);

temp = t-s;
for(i=sp->order; i<sp->n3; i++)
	sp->knot[2][i] = s + ((i-sp->order+1)*temp)/(sp->n3 - sp->order+1);

// print for verification
/*for(i=0;i<sp->order +sp->n1;i++)
	printf("%f, ",sp->knot[0][i]);
printf("\n");
for(i=0;i<sp->order +sp->n2;i++)
	printf("%f, ",sp->knot[1][i]);
printf("\n");
for(i=0;i<sp->order +sp->n3;i++)
	printf("%f, ",sp->knot[2][i]);*/

}

struct TP_2D_DATA * data_2d_initialize(int lenx, int leny, int  * inform) 
                                 
{
/*  Purpose: This function allocates and initializes all the pointers
    in the TP_2D_DATA data structure.  

    Return Value:  A pointer to the newly initialized TP_2D_DATA structure
	               In particular space is allocated so that
				   data.val[i][j] is set up to store the value at (x_i, y_j) 


    Arguments:

        lenx            number of x-values (input)
        leny            number of y-values (input)
        
        
        inform          Error code (output)


    Fatal Errors (which return a NULL pointer):

        inform          Error Condition

          1             lenx  <= 0
          2             leny  <= 0
          4             allocation failure for some member of SP_2D_TP

    */

    struct TP_2D_DATA * dat;
    int i;
    double * temp_ptr1;
    

    /* Initial error testing */
    *inform = 0;
    if(lenx <= 0)
        *inform = 1;
    else if(leny <= 0)
        *inform = 2;
    
    if(*inform > 0) return NULL;



            /*  Allocate space for the SP_2D_TP structure  */
    dat = (struct TP_2D_DATA *) CALLOC(1, sizeof(struct TP_2D_DATA));
    if((void *)dat == NULL){
        *inform = 4;
        return NULL;
    }


            /*  Place numerical values in SP_2D_TP structure  */
    dat->lenx       = lenx;
    dat->leny	    = leny;
    

            /*  Allocate space for the values
                and set the TP_2D_DATA pointers  */
    temp_ptr1   = (double *) CALLOC(lenx*leny, sizeof(double));
    if((void *)temp_ptr1 == NULL){
        *inform = 4;
        return NULL;
    }

	dat->val = (double **) CALLOC(lenx, sizeof(double *));
    if((void *)(dat->val) == NULL){
        *inform = 4;
        return NULL;
    }

    for(i=0; i< lenx; i++)
        (dat->val)[i] = &(temp_ptr1[i*leny]);

            /*  Allocate space for the x and y values
                and set the TP_2D_DATA pointers  */
   
    
    	temp_ptr1   = (double *) CALLOC(lenx+leny+2 , sizeof(double));
    if((void *)temp_ptr1 == NULL){
        *inform = 4;
        return NULL;}

	dat->x = temp_ptr1;
	dat->y = &(temp_ptr1[1+lenx]);
      
     temp_ptr1   = (double *) CALLOC(lenx+leny+2 , sizeof(double));
    if((void *)temp_ptr1 == NULL){
        *inform = 4;
        return NULL;}

	dat->weightx = temp_ptr1;
	dat->weighty = &(temp_ptr1[1+lenx]);

	for(i=0; i< lenx; i++)
		dat->weightx[i] = 1.;

	for(i=0; i< leny; i++)
		dat->weighty[i] = 1.;

	return dat;

}




double ** tp_2d_value(struct SP_2D_TP * sp, double *x, 
				  int lenx, double *y, int leny, int der1, int der2) 
{
/*  Purpose: This function computes the values of a 2d tensor product spline
             or one of its derivatives on a rectangular grid of points. 

    Return Value:  A doubly indexed double array (call it mat) with 
	               mat[i][j] = sp_2d_value(sp, x[i], y[j], der1, der2). 
				   This array is lenx by leny.
				   



    Arguments:

        sp				2d tensor product spline (input)
		x               vector of x-values(input)
		lenx			length of the vector x (input)
        y               vector of y-values (input)
		leny            length of the vector x (input)
        der1			the order of the partial with respect to x (input)
		der2            the order of the partial with respect to y (input)
        
        There are no Fatal Errors.

        

    */
	int i, j, p, q, k, n1, n2, left, mflag;
	int *temp4, *temp5;
    double **temp1, **temp2, **mat, *biatx, sum;

	k = sp->order;
	n1 = sp->n1;
	n2 = sp->n2;
    
	temp4   = (int *) CALLOC(lenx, sizeof(int));
	temp5   = (int *) CALLOC(leny, sizeof(int));
	biatx   = (double *) CALLOC( 
		      (1+sp->order)*(1+sp->order), sizeof(double));

    mat = get_mat(lenx, leny);
	temp1 = get_mat(lenx, k);
	temp2 = get_mat(leny, k);

	for(i=0; i<lenx; i++){
        
        interv(sp->knot[0]-1, n1+k, x[i], &left, &mflag);
		if(mflag == 0){
			bsplvd(sp->knot[0]-1, k, x[i], left, biatx, der1+1);
		    temp4[i] = left;
		    for(j=0; j<k; j++)
				temp1[i][j] = biatx[j+1+k*(der1+1)];
		}
		else{
			for(j=0; j<k; j++)
				temp1[i][j] = 0.; 
		    temp4[i] = k;
		 }
	}
    for(i=0; i<leny; i++){
        
        interv(sp->knot[1]-1, n2+k, y[i], &left, &mflag);
		if(mflag == 0){
			bsplvd(sp->knot[1]-1, k, y[i], left, biatx, der2+1);
		    for(j=0; j<k; j++)
				temp2[i][j] = biatx[j+1+k*(der2+1)];
			temp5[i] = left;
		}
		else{
            for(j=0; j<k; j++)
				temp2[i][j] = 0.; 
			temp5[i] = k;
		}
	}

	for(i=0; i<lenx; i++)
		for(j=0; j<leny; j++){
			sum = 0.;
				for(p=0; p<k; p++){
				for(q=0; q<k; q++)
					sum += sp->coef[temp4[i]-k+p][temp5[j]-k+q]*
					temp1[i][p]*temp2[j][q];}
				mat[i][j] = sum;

		}

free_mat(temp1), free_mat(temp2), free(temp4), free(temp5), free(biatx);
return mat;
}

double sp_1d_value(struct SP_1D * sp, double x, int deriv)
{ 
/*  Purpose: This function computes the values of a 1d spline
             or one of its derivatives at x. 

    Return Value:  The value of the (deriv)-th derivative of the
	               the spline stored in sp at x. 



    Arguments:

        sp				2d tensor product spline (input)
		x               vector of x-values(input)
        y               vector of y-values (input)
		deriv			the order of the derivative of sp at x (input)
		                deriv must be non-negative.
        
        
        There are no Fatal Errors.

        

    */	
	return bvalue(sp->knot -1, sp->coef-1, sp->n, sp->order,
		          x, deriv);
}
	
double sp_2d_value(struct SP_2D_TP * sp, double x,	double y, int der1, int der2)
{
/*  Purpose: This function computes the values of a 2d tensor product spline
             or one of its derivatives at the point (x,y). 

    Return Value:  The spline value
				   



    Arguments:

        sp				2d tensor product spline (input)
		x               x value(input)
        y               y value (input)
		der1			the order of the partial with respect to x (input)
		der2            the order of the partial with respect to y (input)
                        Both der1 and der2 should be nonnegative.
        
        There are no Fatal Errors.

        

 */
	int i, left, mflag;
   double *temp, temp2;

   temp  = (double *) CALLOC(sp->order, sizeof(double));
   interv(sp->knot[0]-1, sp->n1+sp->order, x, &left, &mflag);
   if(mflag != 0)return 0.;
	   
   for(i=0; i<sp->order; i++)
	   temp[i] = bvalue(sp->knot[1] -1, sp->coef[i+left-sp->order]-1, sp->n2, sp->order,
		          y, der2);
   temp2 =  bvalue(sp->knot[0] +left-sp->order-1, temp-1, sp->order, sp->order,
		          x, der1);
   free(temp);
   return temp2;
}

double sp_3d_value(struct SP_3D_TP * sp, double x,	double y, double z,
				   int der1, int der2, int der3)
{
/*  Purpose: This function computes the values of a 2d tensor product spline
             or one of its derivatives at the point (x,y). 

    Return Value:  The spline value
				   



    Arguments:

        sp				2d tensor product spline (input)
		x               x value(input)
        y               y value (input)
		der1			the order of the partial with respect to x (input)
		der2            the order of the partial with respect to y (input)
                        Both der1 and der2 should be nonnegative.
        
        There are no Fatal Errors.

        

 */
	int i, j, leftx, lefty, mflag;
	double *temp, *temp1, temp2;

	temp  = (double *) calloc(sp->order,sizeof(double));
	temp1 = (double *) calloc(sp->order,sizeof(double));
	interv(sp->knot[0]-1, sp->n1+sp->order, x, &leftx, &mflag);
	if(mflag !=0) return 0.;
	interv(sp->knot[1]-1, sp->n2+sp->order, y, &lefty, &mflag);
	if(mflag !=0) return 0.;

	for(i=0;i<sp->order;i++){
		for(j=0; j<sp->order; j++){
			temp[j] = bvalue(sp->knot[2]-1, sp->coef[i+leftx-sp->order][j+lefty-sp->order]-1,
				sp->n3, sp->order, z, der3);
		}
		temp1[i] = bvalue(sp->knot[1]+lefty-sp->order-1, temp-1, sp->order, sp->order, y, der2);
	}

	temp2 =  bvalue(sp->knot[0]+leftx-sp->order-1, temp1-1, sp->order, sp->order,
				  x, der1);
	free(temp); free(temp1);
	return temp2;
}


double ** get_mat(int n1, int n2)
{
/*  Purpose: This function allocates space for a n1 by n2 doubly
             indexed array of doubles.

    Return Value:  n1 by n2 matrix implemented as pointer to pointer 
	               to double.  The indexing starts at [0][0] to 
				   [n1-1][n2-1].  And the storage is allocated contiguously.
	               
				   

    Arguments:

        n1				row index (input)
		n2              column index (input)
        
        
        There are no Fatal Errors.

	Related Functions:  free_mat

 */
	int i;
	double ** mat, * temp_ptr;

            /*  Allocate space for the array  */
    temp_ptr   = (double *) CALLOC(n1*n2, sizeof(double));
    if((void *)temp_ptr == NULL){
       /* *inform = 4; */
        return NULL;
    }

	mat = (double **) CALLOC(n1, sizeof(double *));
    if((void *)(mat) == NULL){
        /**inform = 4;*/
        return NULL;
    }

    for(i=0; i< n1; i++)
        mat[i] = &(temp_ptr[i*n2]);

	return mat;
}

double *** get_mat_3d(int n1, int n2, int n3)
{
/*  Purpose: This function allocates space for a n1 by n2 doubly
             indexed array of doubles.

    Return Value:  n1 by n2 matrix implemented as pointer to pointer 
	               to double.  The indexing starts at [0][0] to 
				   [n1-1][n2-1].  And the storage is allocated contiguously.
	               
				   

    Arguments:

        n1				row index (input)
		n2              column index (input)
        
        
        There are no Fatal Errors.

	Related Functions:  free_mat

 */
	int i;
	double *** mat;

	mat = (double ***) CALLOC(n1, sizeof(double **));
    if((void *)(mat) == NULL){
        /**inform = 4;*/
        return NULL;
    }

    for(i=0; i< n1; i++)
        mat[i] = get_mat(n2,n3);

	return mat;
}


void free_mat(double ** mat)
{
/*  Purpose: This function frees space allocated by get_mat.

    Return Value:  void
	               
				   

    Arguments:

        mat				a matrix allocated by get_mat (input)
        
        
        There are no Fatal Errors.

	Related Functions:  get_mat

 */
	free(*mat);
	free(mat);
}

void free_mat_3d(double *** mat)
{
/*  Purpose: This function frees space allocated by get_mat.

    Return Value:  void
	               
				   

    Arguments:

        mat				a matrix allocated by get_mat (input)
        
        
        There are no Fatal Errors.

	Related Functions:  get_mat

 */
		free(**mat);
		free(*mat);
		free(mat);
}

void sp_2d_graph(struct SP_2D_TP * sp, char * file_name,
				 double a, double b, double c, double d,
					  int der1, int der2, int num_graph )
{
/*  Purpose: This function computes the values of a 2d tensor product spline
             (or one of its derivatives) on a num_graph by num_graph grid and
			 writes the results to a file in tab delimited form for easy import
			 into Excel for example. This allows one to obtain a graph of the
			 spline over the region [a, b]X[c, d].  We assume that a < b and 
			 c < d.

    Return Value:  void
				   

    Arguments:

        sp				2d tensor product spline (input)
		file_name       the name of the file to which the data
		                will be written (input)
        a               Left end point of the x-interval (input)
		b               Right end point of the x-interval (input)
		c               Left end point of the y-interval (input)
		d               Right end point of the y-interval (input) 
		der1			the order of the partial with respect to x (input)
		der2            the order of the partial with respect to y (input)
                        Both der1 and der2 should be nonnegative.
        num_graph       The number of points in the x and y direction (input)
        
        There are no Fatal Errors.

 */
int i, j;
double *x, *y, **values, temp;
FILE *fp;

fp = fopen(file_name, "w");

     x = (double *) CALLOC(num_graph, sizeof(double));
	 y = (double *) CALLOC(num_graph, sizeof(double));
	 
temp = (b-a)/(num_graph-1.);
    for(i=0; i<num_graph; i++) x[i] = a + i*temp;

temp = (d-c)/(num_graph-1.);
	for(i=0; i<num_graph; i++) y[i] = c + i*temp;

 values = tp_2d_value(sp, x, num_graph, y, num_graph, der1, der2);
		    
 for(i=0; i<num_graph; i++){
	 for(j=0; j<num_graph; j++)
		 fprintf(fp, "%7.5f \t", values[i][j]);
	 fprintf(fp, "\n");
 }

free_mat(values), free(x), free(y);  
fclose(fp);
}

double test_eval(struct SP_2D_TP * sp, int der1 , int der2 )
{
/* This is an internal test routine to verify that there are no
indexing or other anomolies in the evaluation routines.  This routine
evaluates the derivative of a spline two different ways, takes the 
square of the difference, and adds up these deviations.  If there are
no errors test_eval should always return zero.  If it does not then there
is an error in the evaluation routines.
*/

int i, j, NUM_GRAPH = 13;
double *x, **values, sum;


     x = (double *) CALLOC(NUM_GRAPH, sizeof(double));
	 

 for(i=0; i<NUM_GRAPH; i++) x[i] = i/(NUM_GRAPH -1.);
 values = tp_2d_value(sp, x, NUM_GRAPH, x, NUM_GRAPH, der1, der2);

 sum = 0.;
 for(i=0; i<NUM_GRAPH; i++){
	 for(j=0; j<NUM_GRAPH; j++)
		 sum += pow(values[i][j]-sp_2d_value(sp, x[i], x[j], der1, der2), 2);
 }

free_mat(values), free(x);  
 return sum;

}

struct TP_2D_DATA * tp_2d_read(char *file_name, int option, int *inform)
{
/*  Purpose:    This function reads data from the data file
                specified in file_name.  The first line of the file
                must contain only the two integers lenx and leny.
                

    Return Value:  (TP_2D_DATA *) A pointer to a Least Squares Tensor
							   Product Data set


    Arguments:

        file_name       The file_name of the data set (input)
        option          Option describing the format of the
                        file to be read(input)
                        {
                            option =  0

                        Data in TP format
						4 2
                        1.1 1.2 1.3 1.4
						2.1 2.2 2.3 2.4

                         
                        }

        inform          Error code (output)


        Fatal Errors:

        inform          Error Condition

          1             set_io_file failed to open the designated
                        input file

          2             data_2d_initialize failed to allocate the 
                        TP_2D_DATA data structure.

          3             set_io_file failed to close the input file


*/

    struct TP_2D_DATA * dat;
    int i, j, lenx, leny, info;
    double z;
	FILE * FIN;

    FIN  = fopen(file_name, "r");
    if( (int)FIN == 0){
        *inform = 1;
        return NULL;
    }

    fscanf(FIN, "%d %d  \n", &lenx, &leny);

    dat = data_2d_initialize(lenx, leny, &info);
    if(info != 0){
        *inform = 2;
        return NULL;
    }
    if(option == 0){
        /*  Data in TP format
			4 2
            1.1 1.2 1.3 1.4
			2.1 2.2 2.3 2.4
              */
        for(j=0; j<leny; j++)
            for(i=0; i<lenx; i++){
                fscanf(FIN, "%lf ", &z );
                dat->val[i][j] = (double) z;
            }
         for(i=0; i<lenx; i++) dat->x[i] = i/(lenx-1.);
		 for(i=0; i<leny; i++) dat->y[i] = i/(leny-1.);
    }
    else{
		fclose(FIN);
        return NULL;
    }
    fclose(FIN);
    
    return dat;
}

double sp_2d_error(struct SP_2D_TP * sp, struct TP_2D_DATA * dat)
{
/*  Purpose: This function returns the raw sum of squares between
             the spline sp and the data values stored in dat->val 
			 on the grid (dat->x) X (dat->y).


    Return Value:  The sum of squres error (SSE)
				   
    Arguments:

        sp				2d tensor product spline (input)
		dat             2d tensor product data set (input)
        
    
        There are no Fatal Errors.

 */
    int i, j;
	double ** values, sum;
	values = tp_2d_value(sp, dat->x, dat->lenx, dat->y, dat->leny, 0, 0);
	sum = 0.;

	
	for(i=0; i<dat->lenx; i++)
		for(j=0; j<dat->leny; j++)
			sum += pow(values[i][j] - dat->val[i][j], 2);

    free_mat(values);
	return sum;
}





void get_max_min(struct TP_2D_DATA *dat, double * max, double * min)
{
/*  Purpose: This function returns the max and min of the data stored
             in dat.


    Return Value:  void
				   
    Arguments:

        dat				a pointer to the TP_2D_DATA structure (input)
		max             the maximum of the data (output)
		min             the minimum of the data (output)
        
    
        There are no Fatal Errors.

 */
	int i, j;
	double temp;
    
	*max = dat->val[0][0], *min = dat->val[0][0];

	for(i=0; i<dat->lenx; i++)
		for(j=0; j<dat->leny; j++){
			temp = dat->val[i][j];
			if(*max < temp) *max = temp;
			if(*min > temp) *min = temp;
		}
}


double plus(double x, int k)
{
/*  Purpose: This function returns (x)_+^k,  that is, 
             x^k if x>0 or 0 if x<= 0; 


    Return Value:  (x)_+^k
				   
    Arguments:

        x				argument for the "+" function (input)
		k               power for the "+" function (input)
        
    
        There are no Fatal Errors.

 */
	if(x > 0) return pow(x, k);
	return 0.;
}


void  bchslv ( double **w, int nbands, int nrow, double *b )
{/*
c  from  * a practical guide to splines *  by c. de boor    
c  solves the linear system     c*x = b   of order  n r o w  for  x
c  provided  w  contains the cholesky factorization for the banded (sym-
c  metric) positive definite matrix  c  as constructed in the subroutine
c    b c h f a c  (quo vide).
c
c******  i n p u t  ******
c  nrow.....is the order of the matrix  c .
c  nbands.....indicates the bandwidth of  c .
c  w.....contains the cholesky factorization for  c , as output from
c        subroutine bchfac  (quo vide).
c  b.....the vector of length  n r o w  containing the right side.
c
c******  o u t p u t  ******
c  b.....the vector of length  n r o w  containing the solution.
c
c******  m e t h o d  ******
c  with the factorization  c = l*d*l-transpose  available, where  l  is
c  unit lower triangular and  d  is diagonal, the triangular system
c  l*y = b  is solved for  y (forward substitution), y is stored in  b,
c  the vector  d**(-1)*y is computed and stored in  b, then the triang-
c  ular system  l-transpose*x = d**(-1)*y is solved for  x (backsubstit-
c  ution).
*/
      int  j, jmax, n, nbndm1;
      if (nrow > 1)                  goto label_21;
      b[1] = b[1]*w[1][1];
                                        return;

/*
c                 forward substitution. solve l*y = b for y, store in b.
*/
label_21: nbndm1 = nbands - 1;
     
	for(n=1; n<=nrow; n++){
         jmax = min_int(nbndm1,nrow-n);
         if (jmax < 1)               goto label_30;
		 for(j=1; j<= jmax; j++){
            b[j+n] = b[j+n] - w[j+1][n]*b[n];
		}
label_30:     n=n;
		  }
/*
c     backsubstitution. solve l-transp.x = d**(-1)*y  for x, store in b.
*/
			n = nrow;    
label_39:   b[n] = b[n]*w[1][n];  
         jmax = min_int(nbndm1,nrow-n);
         if (jmax < 1)               goto label_40;
			 for(j=1; j<=jmax; j++){
			      b[n] = b[n] - w[j+1][n]*b[j+n];
			 }
label_40:    n = n-1; 
         if (n > 0)                  goto label_39;
return;
}


void bchfac ( double **w , int nbands, int nrow, double *diag )
{
/*
c  from  * a practical guide to splines *  by c. de boor    
constructs cholesky factorization
c                     c  =  l * d * l-transpose
c  with l unit lower triangular and d diagonal, for given matrix c of
c  order  n r o w , in case  c  is (symmetric) positive semidefinite
c  and  b a n d e d , having  n b a n d s  diagonals at and below the
c  main diagonal.
c
c******  i n p u t  ******
c  nrow.....is the order of the matrix  c .
c  nbands.....indicates its bandwidth, i.e.,
c          c(i,j) = 0 for i-j .ge. nbands .
c  w.....workarray of size (nbands,nrow)  containing the  nbands  diago-
c        nals in its rows, with the main diagonal in row  1 . precisely,
c        w(i,j)  contains  c(i+j-1,j), i=1,...,nbands, j=1,...,nrow.
c          for example, the interesting entries of a seven diagonal sym-
c        metric matrix  c  of order  9  would be stored in  w  as
c
c
c
c
c
c
c        all other entries of  w  not identified in this way with an en-
c        try of  c  are never referenced .
c  diag.....is a work array of length  nrow .
c
c******  o u t p u t  ******
c  w.....contains the cholesky factorization  c = l*d*l-transp, with
c        w(1,i) containing  1/d(i,i)
c        and  w(i,j)  containing  l(i-1+j,j), i=2,...,nbands.
c
c******  m e t h o d  ******
c   gauss elimination, adapted to the symmetry and bandedness of  c , is
c   used .
c     near zero pivots are handled in a special way. the diagonal ele-
c  ment c(n,n) = w[1][n] is saved initially in  diag[n], all n. at the n-
c  th elimination step, the current pivot element, viz.  w[1][n], is com-
c  pared with its original value, diag[n]. if, as the result of prior
c  elimination steps, this element has been reduced by about a word
c  length, (i.e., if w[1][n]+diag[n] .le. diag[n]), then the pivot is de-
c  clared to be zero, and the entire n-th row is declared to be linearly
c  dependent on the preceding rows. this has the effect of producing
c   x(n) = 0  when solving  c*x = b  for  x, regardless of  b. justific-
c  ation for this is as follows. in contemplated applications of this
c  program, the given equations are the normal equations for some least-
c  squares approximation problem, diag[n] = c(n,n) gives the norm-square
c  of the n-th basis function, and, at this point,  w[1][n]  contains the
c  norm-square of the error in the least-squares approximation to the n-
c  th basis function by linear combinations of the first n-1 . having
c  w[1][n]+diag[n] .le. diag[n] signifies that the n-th function is lin-
c  early dependent to machine accuracy on the first n-1 functions, there
c  fore can safely be left out from the basis of approximating functions
c     the solution of a linear system
c                       c*x = b
c   is effected by the succession of the following  t w o  calls:
c     call bchfac ( w, nbands, nrow, diag )       , to get factorization
c     call bchslv ( w, nbands, nrow, b )          , to solve for x.
c
*/
      int  i,imax,j,jmax,n;
      double    ratio;
      if (nrow > 1)                  goto label_9;
      if (w[1][1] > 0.) w[1][1] = 1./w[1][1];
                                        return;
/*      store diagonal of  c  in  diag */
label_9: 
	  
		  for(n=1; n<= nrow; n++){
		    diag[n] = w[1][n];
		  }
/* factorization */
    
	for(n=1; n<= nrow; n++){
         if (w[1][n]+diag[n] > diag[n]) goto label_15;
		for(j=1; j<= nbands; j++){
			w[j][n] = 0.;
		}
                                        goto label_20;
label_15:    w[1][n] = 1./w[1][n];
         imax = min_int(nbands-1,nrow - n);
         if (imax < 1)               goto label_20;
         jmax = imax;
		 for(i=1; i<= imax; i++){
            ratio = w[i+1][n]*w[1][n];
			for(j=1; j<= jmax; j++){
				w[j][n+i] = w[j][n+i] - w[j+i][n]*ratio;
			}
            jmax = jmax - 1;
            w[i+1][n] = ratio;
		 }
label_20:    n=n; 
	}
return;
    
}

void  get_L2_1D_spline( struct L2_1D_DATA * data, struct SP_1D * sp )
{
/* l2appr ( t, n, k, q, diag, bcoef )
c  from  * a practical guide to splines *  by c. de boor    
c  to be called in main program  l 2 m a i n .
calls subprograms  bsplvb, bchfac/slv
c
constructs the (weighted discrete) l2-approximation by splines of order
c  k  with knot sequence  t(1), ..., t(n+k)  to given data points
c  ( tau(i), gtau(i) ), i=1,...,ntau. the b-spline coefficients
c  b c o e f   of the approximating spline are determined from the
c  normal equations using cholesky's method.
c
c******  i n p u t  ******
c  t(1), ..., t(n+k)  the knot sequence
c  n.....the dimension of the space of splines of order k with knots t.
c  k.....the order
c
c  w a r n i n g  . . .  the restriction   k .le. kmax (= 20)   is impo-
c        sed by the arbitrary dimension statement for  biatx  below, but
c        is  n o w h e r e   c h e c k e d   for.
c
c******  w o r k  a r r a y s  ******
c  q....a work array of size (at least) k*n. its first  k  rows are used
c       for the  k  lower diagonals of the gramian matrix  c .
c  diag.....a work array of length  n  used in bchfac .
c
c******  i n p u t  via  c o m m o n  /data/  ******
c  ntau.....number of data points
c  (tau(i),gtau(i)), i=1,...,ntau     are the  ntau  data points to be
c        fitted .
c  weight(i), i=1,...,ntau    are the corresponding weights .
c
c******  o u t p u t  ******
c  bcoef(1), ..., bcoef(n)  the b-spline coeffs. of the l2-appr.
c
c******  m e t h o d  ******
c  the b-spline coefficients of the l2-appr. are determined as the sol-
c  ution of the normal equations
c     sum ( (b(i),b(j))*bcoef(j) : j=1,...,n)  = (b(i),g),
c                                               i = 1, ..., n .
c  here,  b(i)  denotes the i-th b-spline,  g  denotes the function to
c  be approximated, and the  i n n e r   p r o d u c t  of two funct-
c  ions  f  and  g  is given by
c      (f,g)  :=  sum ( f(tau(i))*g(tau(i))*weight(i) : i=1,...,ntau) .
c  the arrays  t a u  and  w e i g h t  are given in common block
c   d a t a , as is the array  g t a u  containing the sequence
c  g(tau(i)), i=1,...,ntau.
c  the relevant function values of the b-splines  b(i), i=1,...,n, are
c  supplied by the subprogram  b s p l v b .
c     the coeff.matrix  c , with
c           c(i,j)  :=  (b(i), b(j)), i,j=1,...,n,
c  of the normal equations is symmetric and (2*k-1)-banded, therefore
c  can be specified by giving its k bands at or below the diagonal. for
c  i=1,...,n,  we store
c   (b(i),b(j))  =  c(i,j)  in  q(i-j+1,j), j=i,...,min_int(i+k-1,n)
c  and the right side
c   (b(i), g )  in  bcoef(i) .
c  since b-spline values are most efficiently generated by finding sim-
c  ultaneously the value of  e v e r y  nonzero b-spline at one point,
c  the entries of  c  (i.e., of  q ), are generated by computing, for
c  each ll, all the terms involving  tau(ll)  simultaneously and adding
c  them to all relevant entries.
*/
      
	  const int kmax = 20;
	  int k, n, i, j, jj, left, leftmk, ll, mm, ntau;
	  double * bcoef, *diag, ** q, *t, biatx[21], dw;
	  double *weight, *tau, *gtau;
		        
      n = sp->n;
      k = sp->order;
	  t = sp->knot -1;
	  ntau = data->nobs;
	  tau = data->xdata -1;
	  gtau = data->ydata -1;
	  weight = data->weight -1;
	  q = get_mat(k+1, n+1);
	  diag =  (double *) CALLOC(n+1, sizeof(double));
      bcoef = (double *) CALLOC(n+1, sizeof(double));
      
      left = k;
      leftmk = 0;
	  for(ll=1; ll<= ntau; ll++){
      
   /*   locate  l e f t  s.t. tau(ll) in (t(left),t(left+1)) */
label_10:   if (left == n)            goto label_15;
            if (tau[ll] < t[left+1])  goto label_15;
            left = left+1;
            leftmk = leftmk + 1;
			goto label_10;
label_15:   bsplvb ( t, k, 1, tau[ll], left, biatx );
/*
c        biatx(mm) contains the value of b(left-k+mm) at tau(ll).
c        hence, with  dw := biatx(mm)*weight(ll), the number dw*gtau(ll)
c        is a summand in the inner product
c           (b(left-k+mm), g)  which goes into  bcoef(left-k+mm)
c        and the number biatx(jj)*dw is a summand in the inner product
c           (b(left-k+jj), b(left-k+mm)), into  q(jj-mm+1,left-k+mm)
c        since  (left-k+jj) - (left-k+mm) + 1  =  jj - mm + 1 .
*/
         
		for(mm=1; mm<=k; mm++){
            dw = biatx[mm]*weight[ll];
            j = leftmk + mm;
            bcoef[j] = dw*gtau[ll] + bcoef[j];
            i = 1;
         
				for(jj=mm; jj<=k; jj++){
                  q[i][j] = biatx[jj]*dw + q[i][j];
                  i = i + 1;
				}
			 }
	  }
/*
c
c             construct cholesky factorization for  c  in  q , then use
c             it to solve the normal equations
c                    c*x  =  bcoef
c             for  x , and store  x  in  bcoef .
*/
      bchfac ( q, k, n, diag );
      bchslv ( q, k, n, bcoef );
	  for(i=0; i<sp->n; i++) 
		  sp->coef[i] = bcoef[i+1];
	  free_mat(q), free(diag), free(bcoef);
return;
}


int min_int(int i, int j)
{int k;
	k = i;
	if(k > j) k = j;
	return k;
}
	
struct L2_1D_DATA * data_1d_initialize(int nobs, int  * inform) 
                                 
{
/*  Purpose: This function allocates and initializes all the pointers
    in the L2_1D_DATA data structure.  As well as defaulting the weights
	to the value 1.

    Return Value:  A pointer to the newly initialized TP_2D_DATA structure


    Arguments:

        nobs            number of x-values (input)
        
        
        
        inform          Error code (output)


    Fatal Errors (which return a NULL pointer):

        inform          Error Condition

          1             nobs  <= 0
          
          4             allocation failure for some member of L2_1D_DATA

    */

    struct L2_1D_DATA * dat;
    int i;
	double * temp_ptr1;
    

    /* Initial error testing */
    *inform = 0;
    if(nobs <= 0)
        *inform = 1;
    
    
    if(*inform > 0) return NULL;



            /*  Allocate space for the L2_1D_DATA structure  */
    dat = (struct L2_1D_DATA *) CALLOC(1, sizeof(struct L2_1D_DATA));
    if((void *)dat == NULL){
        *inform = 4;
        return NULL;
    }


            /*  Place numerical values in L2_1D_DATA structure  */
    dat->nobs       = nobs;
        

            /*  Allocate space for the values
                and set the L2_1D_DATA pointers  */
    temp_ptr1   = (double *) CALLOC(3*nobs, sizeof(double));
    if((void *)temp_ptr1 == NULL){
        *inform = 4;
        return NULL;
    }

	dat->weight = temp_ptr1;
	dat->xdata  = temp_ptr1 + nobs;
	dat->ydata  = temp_ptr1 + 2*nobs;

	for(i=0; i< nobs; i++)
		dat->weight[i] = 1.;
	
	return dat;

}

struct SP_1D * spline_1d_copy(struct SP_1D * sp)
{
	int i, inform;
	struct SP_1D * copy_sp;

	copy_sp = sp_1d_initialize(sp->order, sp->n, &inform);

	for(i=0; i< sp->order + sp->n ; i++)
		copy_sp->knot[i] = sp->knot[i];
	for(i=0; i< sp->n ; i++)
		copy_sp->coef[i] = sp->coef[i];
	return copy_sp;
}


void get_L2_2d_spline( struct TP_2D_DATA * data, struct SP_2D_TP * sp)
{
      
	  const int kmax = 20;
	  int k, n, i, j, jj, left, leftmk, ll, mm, ntau, second, nrows;
	  double ** bcoef, *diag, ** q, *t, biatx[21], dw;
	  double *weight, *tau, **gtau, **bnewcoef;
		        
      n = sp->n2;
      k = sp->order;
	  t = sp->knot[1] -1;
	  ntau = data->leny;
	  tau = data->y -1;
	  gtau = data->val;
	  weight = data->weighty -1;
	  q = get_mat(k+1, n+1);
	  bcoef = get_mat(data->lenx +1, sp->n2 +1 );
	  nrows = data->lenx;
	  diag =  (double *) CALLOC(n+1, sizeof(double));
      second = 0;
start:     
      left = k;
      leftmk = 0;
	  for(ll=1; ll<= ntau; ll++){
      
   /*   locate  l e f t  s.t. tau(ll) in (t(left),t(left+1)) */
label_10:   if (left == n)            goto label_15;
            if (tau[ll] < t[left+1])  goto label_15;
            left = left+1;
            leftmk = leftmk + 1;
			goto label_10;
label_15:   bsplvb ( t, k, 1, tau[ll], left, biatx );
/*
c        biatx(mm) contains the value of b(left-k+mm) at tau(ll).
c        hence, with  dw := biatx(mm)*weight(ll), the number dw*gtau(ll)
c        is a summand in the inner product
c           (b(left-k+mm), g)  which goes into  bcoef(left-k+mm)
c        and the number biatx(jj)*dw is a summand in the inner product
c           (b(left-k+jj), b(left-k+mm)), into  q(jj-mm+1,left-k+mm)
c        since  (left-k+jj) - (left-k+mm) + 1  =  jj - mm + 1 .
*/
         
		for(mm=1; mm<=k; mm++){
            dw = biatx[mm]*weight[ll];
            j = leftmk + mm;
			for(jj=1; jj<= nrows; jj++)
				bcoef[jj][j] = dw*gtau[jj-1][ll-1] + bcoef[jj][j];
            i = 1;
         
				for(jj=mm; jj<=k; jj++){
                  q[i][j] = biatx[jj]*dw + q[i][j];
                  i = i + 1;
				}
			 }
	  }
/*
c
c             construct cholesky factorization for  c  in  q , then use
c             it to solve the normal equations
c                    c*x  =  bcoef
c             for  x , and store  x  in  bcoef .
*/
      bchfac ( q, k, n, diag );
	  for(jj=1; jj<= nrows; jj++)
		  bchslv ( q, k, n, bcoef[jj] );
	  
	  
	  if(second == 0){
      bnewcoef = get_mat(sp->n2 , nrows );
	  for(jj=1; jj<=  nrows; jj++)
		  for(i=1; i<= sp->n2; i++)
			  bnewcoef[i-1][jj-1] = bcoef[jj][i];}
	  	  
	  free_mat(q), free(diag);

	  if(second == 0) {
      free_mat(bcoef);
	  n = sp->n1;
      k = sp->order;
	  t = sp->knot[0] -1;
	  ntau = data->lenx;
	  tau = data->x -1;
	  
	  weight = data->weightx -1;
      bcoef = get_mat(sp->n2 +1, sp->n1 +1 );
	  nrows = sp->n2;
	  q = get_mat(k+1, n+1);
	  diag =  (double *) CALLOC(n+1, sizeof(double));
      gtau = bnewcoef;
	  second = 1;
	  goto start;}
		  /*clean up*/
		  for(i=0; i<sp->n1; i++)
			  for(j=0; j<sp->n2; j++)
				  sp->coef[i][j] = bcoef[j+1][i+1];
      
     free_mat(bcoef), free_mat(gtau);
return;
}

int b_center(double *v1, double * v2, double *v3, double * pt, double * u)
{   double det;
    int i;
	det = (v2[0] - v1[0])*(v3[1] - v1[1]) - (v3[0] - v1[0])*(v2[1] - v1[1]);
	if (det == 0){
		u[0] = -1.e30, u[3] = -1.e30, u[2] = -1.e30;
		return -1;
	}
	u[0] = (v2[0] - pt[0])*(v3[1] - pt[1]) - (v3[0] - pt[0])*(v2[1] - pt[1]);
	u[1] = (pt[0] - v1[0])*(v3[1] - v1[1]) - (v3[0] - v1[0])*(pt[1] - v1[1]);
	u[2] = (v2[0] - v1[0])*(pt[1] - v1[1]) - (pt[0] - v1[0])*(v2[1] - v1[1]);
	for (i=0; i<3; i++) 
		u[i] = u[i]/det;
	
	if(u[0] >= 0 && u[1] >= 0 && u[2] >= 0) return 0;
    else 
		return -1;
}

struct PWL_2D * pwl_2d_initialize(int n1, int n2, int  * inform)
{
/*  Purpose: This function allocates and initializes all the pointers
    in the PWL_2D data structure.  

    Return Value:  A pointer to the newly initialized PWL_2D structure


    Arguments:

        n1              Domain dimension (input)
        n2              Range dimension (input)
        
        
        inform          Error code (output)


    Fatal Errors (which return a NULL pointer):

        inform          Error Condition

          2             n1        <= 0
          3             n2        <= 0
          4             allocation failure for some member of PWL_2D

    */

    struct PWL_2D * sp;
    double * temp_ptr1;
    

    /* Initial error testing */
    *inform = 0;
    if(n1 <= 0)
        *inform = 2;
    else if(n2 <= 0)
        *inform = 3;
    if(*inform > 0) return NULL;



            /*  Allocate space for the SP_2D_TP structure  */
    sp = (struct PWL_2D *) CALLOC(1, sizeof(struct PWL_2D));
    if((void *)sp == NULL){
        *inform = 4;
        return NULL;
    }


            /*  Place numerical values in SP_2D_TP structure  */

    sp->n1		    = n1;
    sp->n2		    = n2;
    

	sp->coef = get_mat(n1, n2);
	if((void *)(sp->coef) == NULL){
        *inform = 4;
        return NULL;
    }

	
	
	temp_ptr1   = (double *) CALLOC(n1 , sizeof(double));
    if((void *)temp_ptr1 == NULL){
        *inform = 4;
        return NULL;}

    sp->knot_x = temp_ptr1;

	temp_ptr1   = (double *) CALLOC(n2 , sizeof(double));
    if((void *)temp_ptr1 == NULL){
        *inform = 4;
        return NULL;}

    sp->knot_y = temp_ptr1;
	
	return sp;

}

void free_PWL(struct PWL_2D *pwl)
{
	free_mat(pwl->coef);
	free(pwl->knot_x), free(pwl->knot_y);
	free(pwl);
}

double pwl_2d_value(struct PWL_2D * sp, double x, double y)
{
/*  Purpose: This function computes the values of a 2d piecewise linear 
             spline (x,y). 

    Return Value:  The spline value
				   



    Arguments:

        sp				2d piecewise linear spline (input)
		x               x value(input)
        y               y value (input)
		
        
        There are no Fatal Errors.

        

 */
int left_x, left_y, mflag;
double v1[2], v2[2], v3[2], pt[2], u[3], temp;

	if(sp->knot_x[0] > x || sp->knot_x[sp->n1 -1] < x ||
       sp->knot_y[0] > y || sp->knot_y[sp->n2 -1] < y) return 0.;

    interv ( sp->knot_x -1, sp->n1, x, &left_x,  &mflag );
    interv ( sp->knot_y -1, sp->n2, y, &left_y,  &mflag );
    left_x = left_x -1, left_y = left_y-1;

	v1[0] = sp->knot_x[left_x+1], v1[1] = sp->knot_y[left_y];
    v2[0] = sp->knot_x[left_x],   v2[1] = sp->knot_y[left_y+1]; 
	v3[0] = sp->knot_x[left_x],   v3[1] = sp->knot_y[left_y];
	pt[0] = x, pt[1] = y;

	if(b_center(v1,v2,v3,pt, u) >=0) {
		temp = u[0]*(sp->coef[left_x+1][left_y]) +
			   u[1]*(sp->coef[left_x][left_y+1]) +
			   u[2]*(sp->coef[left_x][left_y]);
		return temp;
	}
    v3[0] = sp->knot_x[left_x+1],   v3[1] = sp->knot_y[left_y+1];
    	if(b_center(v1,v2,v3,pt, u) >=0) {
		temp = u[0]*(sp->coef[left_x+1][left_y]) +
			   u[1]*(sp->coef[left_x][left_y+1]) +
			   u[2]*(sp->coef[left_x+1][left_y+1]);
		return temp;
	}
		return -1.e30;
}

struct SP_2D_TP * sp_2d_controlled_surf(struct TP_2D_DATA * data, int num_pts_per_int,
								int out_x, int out_y, double *** valuep, 
								int * inform)
{
/*  Purpose: This function computes a 2d tensor product spline fit
             as well as the values of the fit  
             on an out_x by out_y (equally spaced grid).  

    Return Value:  A pointer to the TP spline representation.
				   



    Arguments:

        data			A TP least squares data set whose values 
		                are to be INTERPOLATED.  (input)
		num_pts_per_int The number of points to be added to control the fit.
		                (input)
		out_x           Number of x-grid values for output (input)
        out_y           Number of y-grid values for output (input)
		valuep			A pointer to the 2D array of values (output)
		
		inform          Error code (output)


    Fatal Errors (which return a NULL pointer):
        
		inform          Error Condition

          1             data->lenx or data->leny  <= 1
		  2				num_pts_per_int <= 0
		  3				out_x <=1 or out_y <= 1
          
          4             allocation failure
      
        

 */

int i, j,  dim_x, dim_y, info;
double *x, *y, **value;
struct PWL_2D * pwl;
struct TP_2D_DATA * data_new;
struct SP_2D_TP *sp;

*inform = 0;
if (data->lenx <= 1 || data->leny <= 1){
	*inform = 1;
	return NULL;}

if (num_pts_per_int <= 0){
	*inform = 2;
	return NULL;}
dim_x = 4 + 2*(data->lenx -2);
dim_y = 4 + 2*(data->leny -2);

sp = sp_2d_initialize(4, dim_x, dim_y, &info);
if (info != 0){
	*inform = 4;
	return NULL;}   

/*  Fill out knot sequences */
for(i=0; i<4; i++){
	sp->knot[0][i] = data->x[0], sp->knot[1][i] = data->y[0];
    sp->knot[0][dim_x+i] = data->x[data->lenx -1];
	sp->knot[1][dim_y+i] = data->y[data->leny -1];
}

j = 4;
for(i = 1; i< data->lenx -1; i++){
	sp->knot[0][j] = data->x[i];
	j = j+1;
	sp->knot[0][j] = data->x[i];
	j=j+1;
}

j = 4;
for(i = 1; i< data->leny -1; i++){
	sp->knot[1][j] = data->y[i];
	j = j+1;
	sp->knot[1][j] = data->y[i];
	j=j+1;
}

/*  Initialize and fill out PWL function */
pwl = pwl_2d_initialize(data->lenx, data->leny, &info);
if (info != 0){
	*inform = 4;
	return NULL;} 

for(i=0; i<data->lenx; i++)
	for(j=0; j<data->leny; j++)
		pwl->coef[i][j] = data->val[i][j];

for(i=0; i<data->lenx; i++)
	pwl->knot_x[i] = data->x[i];

for(i=0; i<data->leny; i++)
	pwl->knot_y[i] = data->y[i];


data_new = data_2d_initialize((data->lenx-1)*num_pts_per_int+1, 
							  (data->leny-1)*num_pts_per_int+1, &info) ;

if (info != 0){
	*inform = 4;
	return NULL;} 

/* Fill out data in data_new */

data_new->x[data_new->lenx-1] = data->x[data->lenx-1];
for(i=0; i<data->lenx-1; i++)
	for(j=0; j<num_pts_per_int; j++)
		data_new->x[num_pts_per_int*i +j] = 
			data->x[i] + j*(data->x[i+1] - data->x[i])/(num_pts_per_int);

data_new->y[data_new->leny-1] = data->y[data->leny-1];
for(i=0; i<data->leny-1; i++)
	for(j=0; j<num_pts_per_int; j++)
		data_new->y[num_pts_per_int*i +j] = 
			data->y[i] + j*(data->y[i+1] - data->y[i])/(num_pts_per_int);


for(i=0; i<data_new->lenx; i++)
	for(j=0; j<data_new->leny; j++)
		data_new->val[i][j] = pwl_2d_value(pwl, data_new->x[i], data_new->y[j]);

for(i=0; i<data_new->lenx; i++)
	data_new->weightx[i] = 1.e-4;

for(i=0; i<data_new->leny; i++)
	data_new->weighty[i] = 1.e-4;

for(i=0; i<data->lenx; i++)
	data_new->weightx[i*num_pts_per_int] = 1.;

for(i=0; i<data->leny; i++)
	data_new->weighty[i*num_pts_per_int] = 1.;

/* compute least squares fit */
get_L2_2d_spline(data_new, sp);

/* get values on a grid */
if (out_x <= 1 || out_y <= 1){
	*inform = 3;
	return NULL;}


x   = (double *) CALLOC(out_x, sizeof(double));
y   = (double *) CALLOC(out_y, sizeof(double));
	
for(i=0; i<out_x; i++)
	x[i] = data->x[0] + i*(data->x[data->lenx-1] - data->x[0])/(out_x-1.);
for(i=0; i<out_y; i++)
	y[i] = data->y[0] + i*(data->y[data->leny-1] - data->y[0])/(out_y-1.);


value = tp_2d_value(sp, x, out_x, y, out_y, 0, 0); 
*valuep = value;

free(x), free(y); 
free(pwl);				 
return sp;
}

struct SP_1D * sp_1d_controlled_curve(struct L2_1D_DATA * data, 
						int num_pts_per_int, int * inform)
{
/*  Purpose: This function computes a 1d spline fit.
             

    Return Value:  A pointer to the spline representation.
				   



    Arguments:

        data			A L2_1D_DATA least squares data set whose values 
		                are to be INTERPOLATED.  (input)
		num_pts_per_int The number of points to be added to control the fit.
		                (input)
		
		
		inform          Error code (output)


    Fatal Errors (which return a NULL pointer):
        
		inform          Error Condition

          1             data->lenx 
		  2				num_pts_per_int <= 0
		  
          
          4             allocation failure
      
        

 */

int i, j,  dim_x, dim_pwl, info;
struct L2_1D_DATA * data_new;
struct SP_1D *sp, *pwl;

*inform = 0;
if (data->nobs <= 1 ){
	*inform = 1;
	return NULL;}

if (num_pts_per_int <= 0){
	*inform = 2;
	return NULL;}
dim_x = 4 + 2*(data->nobs -2);
dim_pwl = data->nobs;

sp = sp_1d_initialize(4, dim_x, &info);
if (info != 0){
	*inform = 4;
	return NULL;}   

/*  Fill out knot sequences */
for(i=0; i<4; i++){
	sp->knot[i] = data->xdata[0];
    sp->knot[dim_x+i] = data->xdata[data->nobs -1];
}

j = 4;
for(i = 1; i< data->nobs -1; i++){
	sp->knot[j] = data->xdata[i];
	j = j+1;
	sp->knot[j] = data->xdata[i];
	j=j+1;
}



/*  Initialize and fill out PWL function */
pwl = sp_1d_initialize(2, data->nobs, &info);
if (info != 0){
	*inform = 4;
	return NULL;} 

for(i=0; i<2; i++){
	pwl->knot[i] = data->xdata[0];
    pwl->knot[dim_pwl+i] = data->xdata[data->nobs -1];
}


for(i = 1; i< data->nobs -1; i++){
	pwl->knot[1+i] = data->xdata[i];
	
}

for(i = 0; i< data->nobs ; i++){
	pwl->coef[i] = data->ydata[i];
	
}
/* PWL is now defined */

data_new = data_1d_initialize((data->nobs-1)*num_pts_per_int+1, &info);

							 
if (info != 0){
	*inform = 4;
	return NULL;} 

/* Fill out data in data_new */

data_new->xdata[data_new->nobs-1] = data->xdata[data->nobs-1];
for(i=0; i<data->nobs-1; i++)
	for(j=0; j<num_pts_per_int; j++)
		data_new->xdata[num_pts_per_int*i +j] = 
			data->xdata[i] + j*(data->xdata[i+1] - data->xdata[i])/(num_pts_per_int);



for(i=0; i<data_new->nobs; i++)
		data_new->ydata[i] = sp_1d_value(pwl, data_new->xdata[i], 0);

for(i=0; i<data_new->nobs; i++)
	data_new->weight[i] = 1.e-6;

for(i=0; i<data->nobs; i++)
	data_new->weight[i*num_pts_per_int] = 1.;



/* compute least squares fit */
get_L2_1D_spline(data_new, sp);


return sp;
}


void free_L2_1D_DATA(struct L2_1D_DATA *p)
{
	free(p->weight);
	/*free(p->ydata); 
	free(p->xdata); 
	
	these statements are unnecessary 
	due to the allocation mechanisim*/
	free(p);
}

void free_SP_1D(struct SP_1D * s)
{
	free(s->coef);
	free(s->knot);
	free(s);
}

void free_SP_2D_TP(struct SP_2D_TP *s)
{
	free_mat(s->coef);
	free_mat(s->knot);
	free(s);
}

void free_SP_3D_TP(struct SP_3D_TP *s)
{
	free_mat(s->knot);
	free_mat_3d(s->coef);
	free(s);
}
